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A lattice fermion model is formulated in Fock space using the Jordan- 
Wigner representation for the fermion creation and annihilation opera- 
tors. The resulting path integral is a sum over configurations of lattice 
site occupation numbers n(x, t) — 0,1 which may be viewed as bosonic 



Ising-like variables. However, as a remnant of Fermi statistics a nonlocal 
sign factor arises for each configuration. When this factor is included in 
measured observables the bosonic occupation numbers interact locally, 
and one can use efficient cluster algorithms to update the bosonized 
rS ' variables. 
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Numerical simulations of lattice fermions are important both in particle and in 
condensed matter physics. For example, in particle physics one simulates quarks 
to investigate if nonperturbative QCD is the correct theory of strong interactions. 
In condensed matter physics one simulates electrons to test if the Hubbard model 
and its variants describe high T c superconductivity. In both cases it is common to 
integrate out the fermions. The resulting bosonic effective theories have nonlocal 
interactions due to the fermion determinant. Because of the nonlocality the theory 
is not really bosonized, and its numerical simulation is computationally difficult 
and very time consuming. Of course, one could include the fermion determinant in 
measured observables and view the system as a local bosonic theory. In practice, 
however, this does not work because the fermion determinant varies over many 
orders of magnitude thus making a numerical simulation extremely inefficient. 

Jordan and Wigner were the first to realize that fermions can be represented 
by bosonic operators Jl|. Their observation was applied to numerical simulations of 
fermionic systems in 1 + 1 dimensions @. The resulting path integral is a sum over 
configurations of occupation numbers n(x, t) =0,1 which may be viewed as bosonic 
Ising-like variables. When generalized to higher dimensions unpleasant minus-signs 
arise which may cause problems in numerical simulations. Still, Duncan || could 
simulate a model of staggered fermions in 2 + 1 dimensions. Later Montvay Q 
refined the method and applied it to Wilson fermions in 3 + 1 dimensions. Here I 
present a new construction which eliminates several constraints present in the con- 
figurations of the previous approaches. As before a nonlocal sign factor arises for 
the Boltzmann weight of each configuration reflecting the Fermi statistics of the 
original theory. Hence, the bosonic theory of occupation numbers is not fully lo- 
cal. However, again one may include this factor in the measured observables and 
view the system without the sign factor as a local bosonized theory. The result- 
ing theory corresponds to a quantum spin system. Of course, also the sign factor 
may fluctuate but, as opposed to the usual fermion determinant, the fluctuations 
are restricted to ±1. Still, for some models this may lead to a minus-sign problem 
of large cancellations between configurations with positive and negative Boltzmann 
weights. In previous approaches to the problem the system of occupation numbers 
was updated using standard local Metropolis or heat bath algorithms. These algo- 
rithms are known to suffer from critical slowing down. Recently, Evertz, Lana and 
Marcu || have developed cluster algorithms for vertex models, which can also be 
applied to quantum spin systems. Unlike the usual Swendsen- Wang- Wolff clusters 
U the Evertz-Lana-Marcu clusters form closed loops. Since our fermion system 
corresponds to a quantum spin system a cluster algorithm can be applied also here. 

To be specific I restrict myself to a simple model which, however, shows the 
characteristic features of more complicated (and physically much more interesting) 
fermionic systems. It is straight forward (although in some cases perhaps tedious) 
to apply the same ideas to more general fermionic models. I consider fermions 
living on the sites of a spatially 2-dimensional L x L lattice with even L and with 



periodic boundary conditions. The characteristic difficulties of fermion updating 
arise because the operators c+ and c x , which create and annihilate fermions at the 
lattice site x = (xi,x 2 ) G Z 2 , anticommute 

(4. C J} = 0, {^, c y } = 0, {c+, c y } = 5 xy . (1) 

The difficulties are due to Fermi statistics, not due to spin. For simplicity I therefore 
consider fermions without spin and with a Hamilton operator 

H = E(4 ^ + c + x+ f x+ - - ctc x+ - - c + x+ f x ), (2) 

x,i 

where i is the unit vector in z-direction. The model is trivial and can be solved 
in momentum space by introducing c+ = \ X^exp^pa;)^, c p = \ J2x ex v{~W x ) c x, 
which implies H = J2 P P 2 Cp c p with pi = 2sin(pi/2). For example, in the grand 
canonical ensemble the expectation value of the occupation number n x = c x c x of a 
site x is given by 

(«.) = ±TV[», eM-W -»N))] = ±Y; 1 + exp()3 1 ( - 2 _ A)) . (3) 

where (3 is the inverse temperature, jjl is the chemical potential and N = J2x n x is 
the fermion number operator. Let us try to obtain the same result from a numerical 
simulation. 

To write the partition function Z as a path integral we first label the lattice 
sites x = (xi,x 2 ) for Xi G {0, 1,...,L — 1} by / G {1, 2, ..., L 2 } in some arbitrary 
order. A convenient choice for numerical applications is e.g. I = 1 + X\ + x 2 L. 
Following Jordan and Wigner |3J the fermion creation and annihilation operators 
are represented as 

c x = (-1) 2 a 1 a 2 ...a l _ 1 a l , c x = (-1) 2 a 1 a 2 ...a l _ 1 a l , (4) 

where a\ are Pauli matrices associated with the lattice point labeled with I and 
a l = \{cr] ± i<?f). The different signs for X\ + x 2 even and odd lattice points are 
not really necessary but they absorb an inconvenient minus-sign in some expressions 
below. I decompose the Hamiltonian into four pieces H = Hi + H 2 + H 3 + if 4 

H\= 22 h x ,i, H 2 = 2_, hx,2, H 3 = 22 h x ,i, H±= 22 ^ 2 > 

x=(2m,n) x=(m,2n) x=(2m+l,n) x=(m,2n+l) 

(5) 
where h xi = ctc x + c + -c^r- — ctc^r- — c + ~.c x . The individual contributions to a 

' a X+l x + l x x + l X+l 

given Hj commute with each other, but two different Hj do not commute. Using 
the Suzuki- Trotter formula one writes for the grand canonical partition function 

Z = Trexp(-(3(H - fiN)) = lim^_ +00 Tr[exp(-e^(F 1 - |jV)) 

x exp(-e/?(tf 2 - !±N)) exp(-e/3(tf 3 - ^N)) exp(-e/3(tf 4 - ^N))] M , (6) 



where e/3 = /3/M is the lattice spacing in the euclidean time direction. Next we insert 
complete sets of Fock states between the factors exp(— e(3{Hj — fiV)). Each site is 
either empty or occupied, i.e. n x has eigenvalue or 1. In the above representation 
by Pauli matrices this corresponds to eigenstates |0) and |1) of af with of |0) = — 10) 
and of |1) = |1). One obtains 



exp(-e/3(h x ,i - -n x - -n x+ -)) = exp(-e/3(l - T )) 



1 n * - i n x+i)) = exp(-e/3(l - - 

/exp(e/?(l -£))A \ 

cosh(e/?)l sinh(e/3)SI 

sinh(e/3)SI cosh(e/3)l 

V exp(-e/3(l - f))l / 



(7) 



where the 4x4 matrix is in the Fock space basis |00), |01), |10), |11) of the sites x 
and x + i. 1 is the identity and SI = of +1 of +2 ...cr^_ 1 is a string of Pauli matrices 
running over consecutive labels between / and m, where / labels x and m labels x + i. 
The operators 1 and SI act on the remaining occupation numbers and are diagonal 
in our basis. The partition function turns into a path integral over occupation 
numbers n(x,t) =0,1 (t labels the time slice) with periodic boundary conditions in 
the euclidean time direction 

Z = U E exp(-5[n])Sign[n]. (8) 

x,t n(x,t)=0,l 

The Boltzmann factor takes the form 

exp(— S[n]) = Yl exj)(—s[n(x,t),n(x + l,t),n(x,t+l),n(x+l,t+l)]) 

x=(2m,n),t=4p 

x Yl exp(—s[n(x,t),n(x + 2,t),n(x,t + i),n(x + 2,t + i)]) 

x=(m,2n),t=4p+l 

x Yl exp(—s[n(x,t),n(x + l,t),n(x,t + l),n(x + l,t + l)]) 

x=(2m+l,n),t=4p+2 

x 'Yl. exp(—s[n(x,t),n(x + 2,t),n(x,t + l),n(x + 2,t + 1)]), 

x=(m,2n+l),t=4p+3 

(9) 

with s[0, 0, 0, 0] = 0, s[l, 1, 1, 1] = 2e/3(l - f ), s[0, 1, 0, 1] = s[l, 0, 1, 0] = e/5(l - 
f) - lncosh(e/?), s[0, 1,1,0] = s[l, 0,0,1] = e/3(l - f ) - lnsinh(e/?). All other 
action values are infinite. Note that the occupation numbers n(x, t) = 0, 1 are 
bosonic variables interacting with each other via the time-like plaquette couplings 
s[n(x, t),n(x + i, t), n(x,t + 1), n(x +i,t + 1)]. This structure is identical with the 
one that occurs in path integral representations of quantum spin systems 0. 

In addition to the Boltzmann factor each configuration is weighted by a sign 
factor which arises from the strings of Pauli matrices of eq. (0) . Just as the Boltzmann 



factor exp(— S[n]) the sign factor Sign[n] is a product of terms siga[n(x,t),n(x + 
i, t), n(x, t + 1), n(x + i,t+ 1)] associated with each plaquette interaction. One has 
sign[0, 0, 0, 0] = sign[l, 1, 1, 1] = sign[0, 1,0, 1] = signfl, 0, 1, 0] = 1. A nontrivial 
sign ±1 may arise only for plaquette interactions of type [0,1,1,0] and [1,0,0,1]. 
To compute the sign factors it is convenient to order the h x { of each Hj in some 
(again arbitrary) way. The factor sign[n(x, t),n(x + i, t),n(x, t + l),n(x + i,t + 1)] 
is the product of eigenvalues ±1 of a^ for the sites with label p between / and m, 
i.e. l + l<p<m — 1. Here I labels the point x and m labels the point x + i. 
If the hyj to which the point with label p belongs is ordered before h x ,i its 4 x 4 
matrix of eq.(^) has already acted on the site labeled by p and we must take its 
occupation number at time t + 1. If, on the other hand, h yi is ordered after h xi 
it has not yet acted on this site and we must take the occupation number at the 
earlier time t. Note that an occupied site gives a factor 1 while an empty site gives 
a factor — 1. Because the individual h X} i of a given Hj commute with each other 
their order does not influence the final result. Moreover, it turns out that the total 
sign factor Sign[n] is independent of the chosen ordering of lattice sites, although the 
contributions from individual plaquettes are in general order-dependent. Hence, any 
reference to the arbitrary order in the Jordan- Wigner representation has disappeared 
from the final expression. This is a very pleasant surprise. The sign factor is 
nonlocal, but it can be computed with an effort proportional to the lattice size if a 
convenient ordering is chosen. More is not needed because the factor is not used in 
the updating process. Since it cannot be interpreted as a probability it is included 
in the measured observables. The system without the sign factor is bosonic and 
interacts locally. In fact, it corresponds to a quantum spin system with Hamiltonian 

H = T, x ,i(i a l a l+i + \ a l a l-H + a l a l + p + ( 2 ~ 2) £* a xi which ma y be viewed as a 
bosonized version of the original fermionic model. This type of bosonization up to a 
nonlocal sign works in any dimension. In ref.|| Ambj0rn and Semenoff proceeded in 
the opposite direction. They fermionized a quantum spin system without the extra 
sign factor and they arrived at 2 + 1-dimensional QED with a Chern-Simons term. 
Luscher also used a Chern-Simons term in a general construction of bosonization in 
2 + 1 dimensions both in the continuum and on the lattice 0. 

The cluster algorithm of Evertz, Lana and Marcu pi, which was originally con- 
structed for vertex models, works very well also for quantum spin systems |10|| . Here 



I describe a variant of the algorithm suitable for updating the fermionic model. The 
algorithm constructs closed loops which are then flipped, i.e. the occupation num- 
bers of points on the loop are changed from to 1 and vice versa. To start a loop 
one first selects a lattice point (x,t) at random. The occupation number n(x,t) 
participates in two plaquette interactions, one at euclidean times before and one at 
euclidean times after t. For n(x,t) = 1 we consider the interaction at the later and 
for n(x,t) = we consider the interaction at the earlier time. The corresponding 
plaquette configuration is characterized by the occupation numbers of four lattice 
points. One of these points will be the next point on the loop. For configurations 
[0, 0, 0, 0] or [1, 1, 1, 1] the next point is with probability p = ~(1 + exp(— ef3)) the 



time-like nearest neighbor of (x, t), and with probability 1— p the next-to- nearest (di- 
agonal) neighbor of (x, t) on the plaquette. For configurations [0,1,0,1] or [1,0,1,0] 
the next point on the loop is with probability pj cosh(e/5) the time-like nearest neigh- 
bor, and with probability 1 — p/ cosh(e(3) the space-like nearest neighbor of (x,t). 
Finally, for configurations [0, 1, 1,0] or [1,0,0, 1] the next point is with probability 
(1 —p)j sinh(e/9) the diagonal neighbor, and with probability 1 — (1 — p)j sinh(e/9) the 
space-like nearest neighbor of (x,t). Once the next point on the loop is determined 
the process is repeated until the loop closes. The above probabilities are arranged 
such that the algorithm obeys detailed balance. The part of the action proportional 
to 1 — j is taken into account by a global Metropolis step. For this purpose each 
loop C is characterized by a winding number W(C), which counts how often the loop 
winds around the lattice in the euclidean time direction. The winding number is 
related to the total occupation number of the loop W(C) = -^ ^2( x ,t)ec(^ n ( x ^ t) ~ 1) 
where 4M is the number of euclidean time slices. The action associates a Boltzmann 
factor exp(— (3(2 — j)W(C)) with each loop. When the loop is flipped its winding 
number changes sign. In the Metropolis step the loop is flipped with probability 
p = min(l, exp((3(4 - fj,)W(C))). 

I have applied the algorithm in a single cluster version for various values of (3 
and \i at fixed lattice spacing e(3 = 1/16. For equilibration 100000 loop clusters have 
been updated, followed by 100000 measurements each separated by 10 loop updates. 
The Monte-Carlo (MC) results are compared to the exact results in table 1; both 
agree within error bars. Note that (n x ) = (n x Sign[n])j,/ (Sign[n])b where b refers to 
the simulated bosonic ensemble. One can see that the minus-sign problem becomes 
more severe when the temperature is lowered or when the chemical potential is 
increased. However, only at (3 — 1, // = 4 there is a real problem, because only 
then (Sign[n])b is consistent with zero. A detailed analysis of the efficiency of the 
algorithm will be presented elsewhere. 

The purpose of the present paper was to describe the bosonization scheme for 
lattice fermions, and to apply it to numerical simulations of a simple model. Of 
course, the idea is to apply the method to models of physical interest. In condensed 
matter physics one can attack the Hubbard model or other models relevant for high 
T c superconductivity. In relativistic lattice field theories one may first study free 
Wilson or staggered fermions. The simplest nontrivial particle physics application 
is perhaps to the purely fermionic Gross-Neveu model. In general, however, the 
fermions are coupled to bosonic fields as well. For example, the Yukawa coupling 
to a scalar field corresponds to a space-time dependent chemical potential and is 
easy to incorporate. In QCD one wants to study the coupling of quarks to SU(2>) 
gauge fields, which should be possible along similar lines. In all cases the minus-sign 
problem might raise its ugly head once one enters a physically interesting regime. 
Only a detailed investigation of the various models of interest can show if this is the 
case or not. Work in some of these directions is in progress. 

I like to thank P. Wiese for a very interesting discussion. 
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s of the numerical simulations. 
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0.976(1) 
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0.25 
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0.4999(7) 
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0.868(2) 


0.1555(5) 


0.15655 


0.50 
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0.779(2) 


0.223(1) 


0.22310 


0.50 
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0.686(2) 


0.306(1) 


0.30493 


0.50 
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8 


32 


0.587(2) 


0.500(3) 


0.50000 


1.00 
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64 


0.747(2) 


0.067(1) 


0.06853 


1.00 


1 


8 


64 


0.322(3) 


0.137(2) 


0.13509 


1.00 


2 


8 


64 


0.052(3) 


0.23(2) 


0.23209 


1.00 


4 


8 


64 


-0.001(3) 


0.5(1.7) 


0.50000 
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